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I. INTRODUCTION 



Antiferromagnetically interacting Heisenberg spins on 
geometrically frustrated lattices have attracted much 
attention of physicists during last yearsii"— A rapidly 
developing direction in this area is the study of the 
properties of such models in the presence of an ex- 
ternal magnetic field. The recent finding that a wide 
class of geometrically frustrated quantum spin anti- 
ferromagnets (including kagome, checkerboard and py- 
rochlore lattices) has quite simple ground states in the 
vicinity of the saturation field - the so-called indepen- 
dent localized-magnon states^i^ - has further stimulated 
studies of the corresponding frustrated quantum anti- 
ferromagnets at high magnetic fields. Among others, 
we mention here the recent papers concerning the de- 
tailed analysis of the low-temperature high-field magne- 
tothermodynamics of a number of one-, two-, and even 
three-dimensional frustrated quantum antiferromagnets 
which support localized-magnon states ^"ii Thus, the 
low-energy degrees of freedom of the quantum Heisen- 
berg antiferromagnet on a kagome lattice in the vicin- 
ity of the saturation field can be mapped onto a gas 
of hard hexagons on a triangular lattice jii^iiiS The lat- 
ter model exhibits a phase transitioni^ that implies a 
phase transition in the spin model at finite (low) tem- 
peratures below (but close to) the saturation field. Sim- 
ilarly, the low-energy degrees of freedom of the checker- 
board antiferromagnet can be mapped onto a gas of hard 
squares on a square lattice^iii with the size of squares 
which corresponds to nearest-neighbor and next-nearest- 
neighbor exclusion and the resulting lattice-gas model 
also exhibits a phase transitions^ that implies the corre- 
sponding peculiarities of the spin model at low tempera- 
tures below the saturation field. Although the performed 
analysis^iiSiii suggests interesting examples of the two- 



dimensional Heisenberg system with a phase transition 
at high magnetic fields and low temperatures, the results 
for the kagome and checkerboard antiferromagnets can- 
not be considered as conclusive examples, since not all of 
the relevant low-energy states are included in the hard- 
core-object lattice-gas descriptioni^iiSrJliil The effect 
of these additional states on the thermodynamic proper- 
ties for both models remains an unresolved problem. In a 
recent paper we have discussed another two-dimensional 
frustrated quantum Heisenberg antiferromagnet - a frus- 
trated bilayer. The low-energy degrees of freedom of 
that model around the saturation field can be mapped 
on a hard-square model (hard squares on a square lat- 
tice corresponding to nearest- neighbor exclusion only), 
see Refs. [l3 and [l3l. Contrary to the kagome and checker- 
board antiferromagnets, for the frustrated bilayer anti- 
ferromagnet the hard-square states completely exhaust 
all low-energy states of the spin model in the vicinity of 
the saturation field and all other low-lying excited states 
are separated by a finite energy gap. Therefore a phase 
transition inherent in the hard-square modeli^ leads to 
a phase transition for the spin model at high magnetic 
fields and low temperatures which cannot be questioned. 
We also note here that in spite of the fact that the 
Mermin- Wagner theorem forbids long-range order for the 
two-dimensional Heisenberg model at any non-zero tem- 
perature at zero field in the presence of an external 
magnetic field it may be indeed present. We will show 
in our paper that emergent discrete degrees of freedom 
may lead to Ising-like antiferromagnetic long-range order 
at low temperatures. 

In the present paper, we extend substantially the stud- 
ies of low-temperature properties for the frustrated bi- 
layer reported in Refs. |l2| and 1 13 and the frustrated two- 
leg ladder reported in Refs.l20l.]lll and[l3lnow taking into 
account within lattice-gas description not only the highly 



degenerate ground-state manifold but also low-lying ex- 
cited states. As a result we arrive at a lattice-gas model 
with finite nearest-neighbor repulsion. For that effective 
model we perform the transfer-matrix calculations (for 
the frustrated ladder) and classical Monte Carlo simu- 
lations (for the frustrated bilayer) to examine the low- 
temperaturc behavior of the quantum spin model for a 
wide region of the magnetic field. Moreover, due to the 
inclusion of the excited states the lattice-gas description 
excellently describes the spin physics up to significantly 
higher temperatures and in a much wider range of the 
magnetic field in comparison with earlier studiesi^iii"— 
The main message of our study is as follows: Geomet- 
rical frustrations may lead to a simple structure of low- 
lying energy levels which can be mapped onto a classical 
lattice- gas model and, as a result, transfer-matrix cal- 
culations or classical Monte Carlo simulations provide 
a very good description of the low-temperature physics 
of the quantum spin model. The most prominent result 
concerns the phase transition in the two-dimensional case 
which "survives" if low-lying excited states are taken into 
account. 

The theoretical investigation of the quantum Heisen- 
berg antiferromagnet on the two-leg ladder and square- 
lattice bilayer has attracted a lot of attention during 
last years. So far the main focus was on ground-state 
properties, see, e.g., Refs. [20l - [26l for the ladder and 
Ref. for the bilayer. In our study we concentrate 
on low-temperature properties of these models. We also 
note that the models under considerations are known as 
models with local conservation laws, see Refs. [20II21II23I - 
[26II28II29I . On the other hand, these models belong 
to a class of systems which support localized-magnon 
states-^iii^i^ 

The paper is organized as follows. In Sec.|ll]wc intro- 
duce the quantum spin models and discuss their symme- 
tries. In Sec. |in]we briefly present a spectroscopic study 
of the spin models. In particular, we focus on a class 
of simple product eigenstates [independent (isolated) lo- 
calized magnons and interacting (overlapping) localized 
magnons] which become the low-energy states under cer- 
tain conditions. In Sec. llVl we map the localized-magnon 
states on the lattice-gas-model states and discuss the de- 
generacy of the localized-magnon states. In Sec. |V] we 
calculate the contribution of the independent localized- 
magnon states to the thermodynamic quantities, whereas 
in Sec. IVII we extend calculation of thermodynamic quan- 
tities taking into account the set of interacting localized- 
magnon states in addition. In these sections we compare 
latticc-gas model results with exact diagonalization data 
for finite spin systems to find the region of validity for 
the lattice-gas-model description. Moreover, we obtain 
the low-temperature thermodynamic quantities of the 
quantum spin models on the basis of (i) transfer-matrix 
calculations for the one-dimensional case and (ii) clas- 
sical Monte Carlo simulations for large two-dimensional 
latticc-gas systems. A brief summary is given in the last 
section (Sec. |VlT|. Some lengthy formulas for one- and 
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FIG. 1: (Color online) Lattices considered in this paper: (a) 
the frustrated two-leg ladder and (b) the frustrated bilayer. 
The vertical bonds have the strength J2 > whereas all other 
bonds have the strength Ji > 0. 



two-dimensional lattice-gas models are collected in two 
appendices. 



II. MODEL 

We consider the antiferromagnetic Heisenberg model 
of = 2M quantum spins s = 1/2 on the two lattices 
shown in Fig. [T] The Hamiltonian of the model reads 

= ^ JpgSp • s, - /i^^ (2.1) 

Here the sum runs over the bonds which connect the 
neighboring sites on the spin lattice shown in Fig. [TJ 
Jpg > are the antiferromagnetic exchange constants be- 
tween the sites p and q which take two values, namely, J2 
for the vertical bonds and J\ for all other bonds, ft, > 
is the external magnetic field, and 5*^ — is the 

z-component of the total spin. In our study we imply 
periodic boundary conditions. Further we set Ji = 1 if 
not stated otherwise explicitly. 

We introduce an underlying lattice oiM ~ N/2 sitesi^ 
For the frustrated two-leg ladder that is a simple chain 
whereas for the frustrated bilayer it is a square lattice. 
Now it is convenient to denote the spin lattice sites as 
TO,i, where m numbers the vertical bonds, i.e., it runs 
over all sites of the underlying lattice and the index i 
refers either to the lower [i = 1) or to the upper [i = 
2) leg or layer. Note that the underlying lattices are 
bipartite ones, i.e., we can divide it into two sublattices 
A and B, and any two neighboring sites on the lattice 
belong to different sublattices. Introducing further the 
total spin of a vertical bond = s^^i -I- Sm.2 the spin 
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Hamiltonian (12.11) can be rewritten as follows 



H = 



2 ' ' 



- - 1 - htf- 



t„.(2.2) 



Here the first sum runs over all sites of the underlying 
lattice and the second sum runs over all bonds which 
connect the neighboring sites on the underlying lattice. 

The Hamiltonian (|2.2p depends on the total spin of 
each vertical bond t,„, m = 1, . . . , A/" only, and the value 
of the total spin of a vertical bond is a good quantum 
number. As a consequence, the properties of the consid- 
ered models can be studied in much more detail. In par- 
ticular, a large number of eigenstates can be constructed 
exactly, see Refs.llMllOand SecUHl In Sec.|lll]we give 
a precise description of low-energy eigenstates of the spin 
Hamiltonian in a regime when J2/J1 is sufficiently large 
(strong-coupling regime). We again emphasize that our 
focus is the low-temperature thermodynamics of these 
models in the strong-coupling regime and therefore we 
will be interested not only in the energies of low-energy 
eigenstates but also in their degeneracies, see Sec. IIVI 



III. PRODUCT EIGENSTATES 

In this section we briefly summarize some known facts 
on a class of simple product eigenstates of the spin Hamil- 
tonian (|2.ip [or (|2.2p ] which become the low-energy ones 
under certain conditions'^"— For this purpose we may 
consider the subspaces with different values of sepa- 
rately since the Hamiltonian (|2.ip commutes with the 
operator . We may assume at first h = since adding 
of the Zeeman contribution is trivial. Obviously, the fully 
polarized state | t, • • ■ , t) is an cigcnstate of the Hamil- 
tonian (EH]) with = N/2 = TV. The energy of this 
state is Epu = J^Ji + A/'J2/4 or Epu = 27V Ji +7VJ2/4 
for the one- or two-dimensional case, respectively. This 
state is the ground state for high magnetic fields. 

Next we consider eigenstates 



\n) = |0,„,)|0„ 
1 



|0™J|FM)r, 



|Omi) — (I tmi,lUni,2) ~ | imi,ltmi,2)) , (3-1) 

where a subset oi n, 1 < n < Af spin pairs on verti- 
cal bonds mi, 7712, m„ are in a singlet state \0m) 
and the remaining M — n other spin pairs on vertical 
bonds are in a fully polarized triplet state | tt,itt,2) with 
tj = 2 and ~ 1, where the index t labels the verti- 
cal bond carrying the triplet. The triplet bonds form a 
fully polarized ferromagnetic background |FM)r. Obvi- 
ously, these states have a magnetization ^ Af — n. 
Each vertical singlet contributes with — 3J2/4 to the en- 
ergy. The contribution of a polarized triplet at a verti- 
cal bond t is J2/4 4- 74^1, where jt counts the number 
of neighboring triplets of the triplet at certain vertical 
bond t. Hence, for the energy of the state (|3.1|) we get 



E„ = -3nJ2/4 4- (TV - n) J2/4 + Ji J2[ 7t/2. It remains 
to calculate X^t Tti where the sum runs over all A/" — n 
vertical triplet bonds. To get a state of minimal energy 
in a certain sector of we have to minimize 7* . 

Obviously, we get the minimal 74 = 2(A/' — 2n) in 
the one-dimensional case or J2[ 7t 4(7V - 2n) in the 
two-dimensional case in all sectors J\f/2 < < J\f, if 
we have no neighboring singlets (hard-core rule). This 
is a weak constraint, and, consequently, there are many 
states fulfilling this constraint. Note that these states 
belong to the class of so-called independent localized- 
magnon states appearing in many frustrated lattices^"— 
The energy of these localized-magnon states is 



E. 



Im 



E, 



FM 



nei. 



(3.2) 



where ei ~ J2 + 2Ji or ei = J2 + 4Ji for the ladder 
or bilayer case, respectively. The maximal number of 
independent localized magnons is n„iax = M/2. For 
n = 'T-max there are two degenerate localized-magnon 
states where one sublattice {A or B) carries the ver- 
tical singlets and the other one is occupied by vertical 
triplets (so-called "magnon-crystal" states). The energy 
of a magnon-crystal state is i?]^/2 ~ — A/'J2/4. In Ref. [ill 
and m (see also Ref. ^l|) it was shown that the indepen- 
dent localized-magnon states are ground states in the re- 
spective sectors of if J2 > 2Ji for the ladder and 
J2 > 4Ji for the bilayer. In a magnetic field we have the 
energy El^{h) = E^"" - hS'' = E^m - nci - h{N/2 - n). 
Hence all these independent localized-magnon states are 
degenerate at h = hi ~ ei, where hi = ei is the satu- 
ration field. As a result one finds for the ground-state 
magnetization AI{T ~ 0,h,N) the well-known jump to 
saturation a,t h = hi with a preceding wide plateau (see, 
e.g., Refs.[3ll3l20[). illustrated in Fig.[2J where the plateau 
state is a two-fold degenerate magnon-crystal state. 

Product eigenstates with higher energies in the sectors 
Af/2 < < Af—1 are states where some of the n vertical 
singlets are neighbors. These states can be understood 
as non-independent (i.e., interacting) localized-magnon 
states. Supposing that we have pairs of neighboring 
vertical singlets then we get an energy 



(3.3) 



where Ji can be understood as the repulsion energy. For 
large enough J2 > eigenstates with v ~ 1 are the low- 
est excitations above the independent localized-magnon 
ground states for M /2 < < Af — 1. Based on finite- 
size calculations {N = 32) we estimate J| w 3.00Ji for 
the ladder and J| w 4.65 Ji for the bilayer. 

Going to lower magnetization < < Af/2 no in- 
dependent localized-magnon states exist, and the "inter- 
acting" localized-magnon states with n > Af/2 vertical 
singlets can become ground states. A lowest-energy state 
with n = Af/2 + r, r = 1, . . . ,N /2 localized magnons 
is, e.g., a state where Af/2 magnons (singlets on verti- 
cal bonds) occupy one sublattice completely (i.e., they 
arc in the "magnon-crystal" state) and the remaining r 
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FIG. 2: (Color online) Ground-state magnetization curves for 
the frustrated ladder (thin solid) and bilayer (thick broken) in 
the considered strong-coupling regime. We set Ji = 1, J2 = 4 
(ladder) and Ji = 1, J2 = 5 (bilayer), i.e., we have /ii = 6 and 
/i2 = 4 for the ladder and /ii = 9 and h2 — 5 for the bilayer. 



magnons sit on the other sublattice. These states have 
a magnetization = f\f/2 — r and an energy (now 

e;7*=o) 



■+r 



rJ2. 



+r 



(3.4) 



hS' = 



In a magnetic field we have the energy Ej^f /2 
— A/'J2/4 — r J2 — h[N /2 — r). Hence all these interacting 
localized- magnon states are degenerate at h = h2 = J2- 
As a result one finds another jump in the ground-state 
magnetization at /i = /12 with a preceding wide = 
plateaup^ see Fig. where this plateau state is a non- 
degenerate state where all vertical bonds carry a sin- 
glet. Low-lying excited states in the subspaces with 
< < N/2 are constructed by rearranging the verti- 
cal singlets to increase the number of neighboring singlets 
to V > 1 (i.e., the sublattice formerly completely occu- 
pied by singlets on vertical bonds is now incompletely 
occupied by singlets). For these excited states the in- 
crease of energy of the resulting states again is given by 
i/Ji, see Eq. 



IV. 



DEGENERACY OF LOCALIZED-MAGNON 
STATES 



After having illustrated the basic facts on exact prod- 
uct eigenstates of the considered models which become 
ground states and excited low-energy states in all sub- 
spaces with = N/2, . . . , for sufficiently strong ver- 
tical bonds J2 > J|, we will now calculate their degen- 
eracies using a mapping of low-energy degrees of free- 
dom of the quantum spin model (j2.ip on appropriate 
classical lattice-gas models. This kind of mapping was 



(a) 



(b) 



FIG. 3: (Color online) Hard-dimer description of the 
localized-magnon states of the frustrated two-leg ladder, (a) 
Independent localized-magnon states (non-overlapping hard- 
dimer states), (b) Interacting localized-magnon states (over- 
lapping hard-dimer states). 



used for various frustrated lattices hosting independent 
localized-magnon statesj^i^^"— Let us mention here that 
recently it has been found that this kind of mapping 
is also applicable to some flat-band Hubbard modelsi^ 
Note that in all previous papers using such a mapping, 
it was restricted to the independent localized-magnon 
ground states, only. Based on the spectroscopic analy- 
sis given in the previous section, here we overcome this 
restriction and extend the results for the ladder and the 
bilayer presented in Refs. fllMlsl including the interacting 
localized-magnon low-lying excited states. 

For a better understanding of the mapping of the inter- 
acting localized-magnon states that will be discussed in 
the next paragraph we start with a brief illustration of the 
mapping of the independent localized-magnon states, see 
also Refs. [TT| - [l3l As mentioned in the previous section 
for the independent localized-magnon states a hard-core 
rule is valid, i.e., they cannot occupy neighboring sites on 
the underlying lattice (chain or square lattice). Hence, 
the number of possibilities to put n independent localized 
magnons on the two-leg ladder (bilayer) is equivalent to 
the number of possibilities to place n hard dimers (hard 
squares) on a chain (square lattice) oi J\f ~ N/2 sites, cf. 
Fig.|3lja). Denoting the ground-state degeneracy in the n- 
magnon subspacc by gj\f{n), we find for n — . . . ,M /2 
that gj\!{n) = Z\ic{n,M), where Z\^c{n,M) is simply the 
canonical partition function of n hard-core objects (hard 
dimers or hard squares) on a A/'-site lattice (chain or 
square lattice) with periodic boundary conditions im- 
posed. We may call the independent localized-magnon 
states the hard-core states. As mentioned above, the in- 
dependent localized-magnon states are ground states in 
the subspaces with 7V/2 < S"" < TV if J2 > 2 Ji (ladder) 
or J2 > 4Ji (bilayer). Moreover, they are linearly inde- 
pendent and form an orthogonal basis in each subspace 
with N/2 < < N, see Ref. i.e., all these states 
contribute to the partition function of the spin system. 
Due to their huge degeneracy they dominate the low- 
temperature thermodynamics for magnetic fields h near 
hi, see Refs. [uHll and Sec. W\ 

Now we extend the above procedure considering the 
interacting localized-magnon states to calculate the de- 
generacy of the low-lying excited states. In the language 
of hard-core objects the hard-core rule is partially re- 
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laxed and the hard-core objects may partially overlap, 
see for illustration Fig. ^h). In what follows we call 
them overlapping hard-core states. Note that a complete 
overlap (i.e., two hard-core objects on the same site) is 
strictly forbidden, since a corresponding spin state does 
not exist. First we consider the lowest excited states in 
the sectors 7V/2 < S"^ < TV - 1 (i.e., > n > 1). 

These states have two and only two neighboring sin- 
glets, i.e., (i) we have only one pair {v = 1) of neighbor- 
ing vertical singlets and (ii) in the hard-core model two 
(and only two) hard-core objects overlap, cf. Fig. [SJb). 
We denote the degeneracy of the first excited states in 
the n-magnon subspace by x^f{n). Clearly, x^{n) for 
n = 2, . . . ,J\f/2 equals the canonical partition function 
of a system with n — 2 non-overlapping hard-core objects 
and one composite hard-core object built by two over- 
lapping objects. In the one-dimensional case we immedi- 
ately conclude, that xj\f{n) = A/'^hc("- — 2, A/"— 4), where 
Zhci'/n, Al) is the canonical partition function of m hard 
dimers on A^-site chain, however, (instead of periodic) 
with open boundary conditions. In the two-dimensional 
case Xf^[n) = 2MZhc{n — 2, A/" — 8), where Zhc(m-, M) is 
the canonical partition function of m hard squares on a 
(periodic) A^-site square lattice with a "dumbbell void" 
oriented either in horizontal or vertical direction, where 
the two overlapping hard squares are located. Using sim- 
ilar reasonings we may find the degeneracy of higher ex- 
cited states with ly > 1. However, since we are interested 
in thermodynamic properties, we need in fact the energies 
and degeneracies of excited states only in certain combi- 
nations which enter the lattice-gas thermodynamics, see 
Sec. ED 

It is straightforward to determine the ground-state 
degeneracy gj\f(ri) in the remaining subspaces, i.e., for 
= Af/2 - 1, . . . , (that is for n = Af/2 + l,...,Af) 
and the first-excited-state degeneracy Xf^in) in the sub- 
spaces with S"^ = N /2 — 1,...,2 (that is for n ~ 
N /2 -\- 1, . . . ,Af — 2). The ground state is then a state 
with minimal overlap of hard-core objects, e.g., one sub- 
lattice is completely occupied and the other one carries 
n—Af/2 hard-core objects. Then the first excited state is 
a state with one (and only one) empty site in one sublat- 
tice together with a neighboring empty site on the other 
sublattice ("composite hole"). As a result, one has a 
simple "particle-hole" symmetry for hard-core objects: 
gj^fin) = gj\f{N' — n) and xj\f{n) = xj\f{N — n). Note that 
such a particle-hole symmetry is valid also for higher ex- 
cited states. 

We can find the degeneracies gj\f{'n) and xj\f{n) ana- 
lytically in the one-dimensional case knowing the parti- 
tion functions Z\^c{n^M) and Zhc('Ti,A4). These quan- 
tities follow from the grand-canonical partition function 
for the one-dimensional hard-dimer model with periodic 
and open boundary conditions through calculation of 
derivatives with respect to the hard-dimer activity z 
at z = 0, e.g., n!Zhc(n,A/') = d"Spbc(z, A/')/(iz"|^=o- 
The grand-canonical partition function for the one- 
dimensional hard-dimer model can be obtained by the 



transfer-matrix method, see, e.g., Ref. [l^ For periodic 
boundary conditions imposed we have 

Spbc(^,AA) =Af + (4.1) 

with A± = (1 ± a/I + ^z)/2. For open boundary condi- 
tions imposed we have 

Sobc(;2, A/") = (a+ + 2VIa+6+ 4- zh\) \^-^ 

+ {a^_ + 2Via-6- + zbi) A^"^ (4.2) 

with the same A± and a± ~ (1 ± a/I + 4^;) / ■y/2C±, 
h± = ^2z/C±, C± = 1 4z ± VI + 4z. In the two- 
dimensional case, the required partition functions can be 
easily found numerically. For example, for M = N/2 = 
16 we get Zhc(n,7V) = 1,16,88,208,228,128,56,16,2 
for n = 0,1,2,3,4,5,6,7,8 or 27VZhc(n - 2,7V- 8) = 
32, 256, 576, 448, 64 for n = 2,3, 4, 5, 6. 

To check the hard-core predictions we have calculated 
the degeneracies of the ground states and the lowest ex- 
citations as well as the excitation gaps of the s = 1/2 
frustrated two-leg ladder and bilayer by full diagonaliza- 
tion for finite spin systems up to = 32 and various 
sectors of total 5*^. The exact diagonalization data coin- 
cide perfectly with the corresponding data obtained by 
analytical formulas JS?!]), (|331), dSZl), (HUl), and (g^) and 
numerics for Zhc(n,A/') and Zhc(m.,A/{). 

As mentioned in Sec. IIIH at the fields hi = ei all 
independent and at /12 — J2 all interacting localized- 
magnon ground states are degenerate which leads to a 
jump in the magnetization curve, see Fig. [2] These 
degeneracies of localized-magnon states at hi and /12, 
Wi = J2n=o 9^fin) and W2 = Y.n=M /2 9M{n), grow ex- 
ponentially with the system size N . Due to the particle- 
hole symmetry we have Wi = Y^2 = W. The exponential 
growth of W leads to a nonzero ground-state residual en- 
tropy S{T = 0,h,N)/N {\nyV)/N ^ at h ^ hi and 
h = h2. 

To summarize Sees. Illll and llVi we have characterized 
the low-energy eigenstates of the frustrated two-leg lad- 
der and bilayer in the strong-coupling regime calculating 
their energies and degeneracies. In the next sections we 
show how due to a simple structure of these low-energy 
eigenstates their contributions to thermodynamics can be 
obtained with the help of auxiliary lattice-gas models. 

V. LATTICE-GAS MODELS WITH 
NEAREST-NEIGHBOR EXCLUSION 
(HARD-CORE MODELS) 

We want to calculate a partition function Z{T, h, N) of 
the spin system (j2.ip . In a first step we consider in this 
section the contribution of independent localized-magnon 
states (non-overlapping hard-core states) to the partition 
function. Recall that these states are the ground states 
in the subspaces with = A^/2, . . . , A^/4 with energies 
El^{h), n = N/2 - 5^ = 0, 1, ... , A^/4 and degeneracies 
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gjsf{n). If the magnetic field h is around the saturation 
filed hi they will give the dominant contribution at low 
temperatures T. Therefore 



n=0 



£gA.(n)e^", (5.1) 



71=0 

where ^ = ei — h = hi — h. Obviously, the magnetic field 
and the temperature enter the thermodynamic quanti- 
ties within the hard-core description via the combination 
(hi — h)/T, only. Since gj^(n) is the canonical partition 
function Ziicin,Af) of n hard dimers on a chain of J\f 
sites or of n hard squares on a square lattice of A/" sites, 
Eiic{T, n,Af) = J2^1o 5Ar('^)e'*"^"^ is the grand-canonical 
partition function of the corresponding onc-dimcnsional 
hard-dimcr model or two-dimensional hard-square model 
and fj. is the chemical potential of the hard-core objects. 
It is also useful to rewrite Z\jn{T, h, N) in the following 
form: 

ZUT,h,N)^e r 

X E •■• E ^-^-""-i^d"™}), (5.2) 

ni— 0,1 n^— 0,1 

where the factor R{{nm}) takes care about the hard- 
core rule, i.e., it is if the spatial configuration {n„i} 
violates the hard-core rule but it is 1 if the hard-core 
rule is fulfilled. For example, for the one-dimensional 
hard dimers i?({n,„}) = (1 — 71171,2) (1 — 7i2fT-3) . . . (1 — 
n^f-ln^f){l - n^/ni). 

In summary, we arrive at the basic relation for the 
independent localized-magnon state contribution to the 
Helmholtz free energy of the spin system (|2.1I) 



FUT,h,N) E, 



FM 



N 



N 



rin5hc(r,At,A/') 
2 M 



(5.3) 



The entropy S, the specific heat C, the (uniform) mag- 
netization M, and the (uniform) susceptibility x follows 
from (|5.3p according to usual formulas, S(T,h,N) = 
-dF{T,h,N)/dT, C{T,h,N) = TdS{T,h,N)/dT, 
M{T,h,N) = N/2 - n, n = Td\nE{T, fi,J^)/dfi, 
X(T, h, N) = dM{T, h, N) /dh = dn/dn. 

To examine the ordering of hard-core objects (local- 
ized magnons) we consider the average total numbers of 
hard-core objects on the sublattices A and B, Wa and ub- 
Obviously n — Wa + ns, whereas the value of the differ- 
ence \nA — tibI may play a role of the order parameter m. 
Introducing an infinitcsimally small symmetry-breaking 
staggered component into the chemical potential, i.e., 
/i — > /i/i = /i + (5/i on the sublattice A and fi ^ I-Ib = 
/i — (5^ on the sublattice B, jj, = hi — h, dfi = —6h, we 
calculate the staggered magnetization Mst{T, h, Sh, N) = 
Ma-Mb = -nX+n^ ^ -Td\nE{T, ^1A, ^J.B,^f)/^^^A + 
Td\nE{T,fiA,HB,^)/dfiB = Xst{T,h,N)6h, where 



Xst = dMst/dSh is the staggered susceptibility. Decreas- 
ing the temperature a divergence of the staggered suscep- 
tibility in the thermodynamic limit signals a transition to 
an ordered phase, where the symmetry of the occupation 
of both sublattices with localized magnons (hard-core ob- 
jects) can be spontaneously broken. 

Below we discuss bricfiy thermodynamic quantities as 
they follow from the lattice-gas models with nearest- 
neighbor exclusion and compare them with exact diag- 
onalization data for the frustrated quantum Heisenberg 
antiferromagnets on finite lattices. 



A. Frustrated two-leg ladder 

We begin with the case of the frustrated two-leg ladder. 
Using the transfer-matrix result for one-dimensional hard 
dimers [sec Eq. (|4.1[) ]. 



Shc(r,/i,AA) = Af + A^, 



1 , /I ii 
A± = -±a/| + ^, z = eT , hi 



(5.4) 



one can easily find all thermodynamic quantities, S, C, 
M, and x, see Eqs. ^M^ - and jM]) - dMl) in the 
Appendix |Al The main features of the low-temperature 
thermodynamic behavior of the frustrated two-leg lad- 
der in this regime are as follows: (i) the jump in zero- 
temperature magnetization at /i = hi, cf. Fig. [21 is 
smeared out at low but finite nonzero temperatures; 
(ii) the entropy S(T,hi,N)/N remains finite and ap- 
proaches S{T = 0,hi,N)/N = (1/2) In ^ as T 0, 
where ip = {1 + V5)/2 is the golden mean [note that due 
to the particle-hole symmetry explained in Sec. IIVI there 
is the same ground-state residual entropy at ft. = /12, 
S{T = 0,h2,N)/N = (1/2) In (iii) the specific heat 
shows an extra low-temperature maximum if h slightly 
deviates from hi indicating a new low-energy scale set- 
tled by the set of independent localized-magnon states. 
A comprehensive analysis of low-temperature high-field 
thermodynamic quantities S, C, M, and x based on 
the hard-dimcr description (j5.3p . (15.41) can be found in 

Rcfs. [Ulil 

To calculated the staggered susceptibility (not consid- 
ered in previous papersii^i^) we have to consider different 
chemical potentials on the sublattices. Then the grand- 
canonical partition function reads 



„ JV AT 

:^hciT,flA,^J■B,^f) = ^+ +c 




(5.5) 



ZA ^ e T , zb ^ e T ^ ^ hi - hA, f^B = hi - hB- 



Eq. ()5.5p immediately yields the staggered susceptibil- 
ity Xst{T,h,N), see Eqs. (|X5|) . (|A10|) . As expected, 
for the one-dimensional problem there is no divergence 
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at T > 0. However, for h < hi the staggered suscep- 
tibility diverges at T = 0. Precisely at ft. = /ii one 
finds Xst{T,hi,N)/M (l/\/5) T'^ « 0.447214 T-\ 
For h < hi one finds Xst{T,h,N)/N = {N/A) T'^ 
for finite M . In the thermodynamic limit A/" — > oo 
the divergence becomes exponential Xst(T,h, N)/Af = 
(l/2)T-ie(''i-'')/(2'r). Note that this temperature de- 
pendence is identical to that of the Ising chain. Note 
further that, trivially, Txst{T,h,N)/Af goes to zero at 
T = Oifh> hi. 



B. Frustrated bilayer 

We turn to the case of the frustrated bilayer. By con- 
trast to the one-dimensional case discussed in the pre- 
vious section there is no exact analytical solution for 
the corresponding two-dimensional hard-square model. 
Thermodynamic functions for the hard-square model can 
be obtained by direct computations only if Af is small 
enough (see Appendix [B|) . For larger Af we use classical 
Monte Carlo simulations^ (see also Appendix [B|) . 

First we illustrate the validity of the hard-square de- 
scription by comparison with exact diagonalization data 
for finite systems. For that we show in Fig. U the specific 
heat C(T, h, N). According to the upper panel of Fig. |4] 
for T < 0.2 the specific heat depends only on the hard- 
core parameter {h — hi) /T and the hard-core description 
is valid over the entire range of magnetic fields. The tem- 
perature dependence of C for various magnetic fields is 
shown in the lower panel of Fig. |4l It is again obvious 
that the low-temperature behavior is well described by 
the hard-square model. However, the temperature range 
of the validity of the hard-square model becomes smaller 
with increasing hi — h. In particular, the position and the 
height of the characteristic extra low-tcmpcrature maxi- 
mum in C is correctly described only if/ii— ft<0.5. 

Similar as for the one-dimensional case at h = hi, there 
is a ground-state residual entropy S{T = 0,hi,N)/N = 
0.2037... = (l/2)lnK(l). This value follows from the 
hard-square model theory, which predicts for the hard- 
square entropy constant k(1) = 1.50304808 . . .i^ More- 
over, due to the particle- hole symmetry, see Sec. lIVi there 
is the same ground-state residual entropy at h ~ /12, 
S{T = 0,h2,N)/N = {1/2)\iik{1). 

The main peculiarity of the low-temperature thermo- 
dynamics of the frustrated bilayer around /ii is connected 
with an order-disorder phase transition which is inherent 
in the hard-square model at Tc « {hi - h)/1.33iOM^ 
The phase transition has pure geometrical origin: if 
the density of hard squares [controlled by the activity 
z = e'^''!"'')/-^] increases they start to occupy only one 
of two sublattices. The critical value of the activity is 
Zc = 3.7962 .... The universality class is that of the two- 
dimensional Ising model, i.e., we have a logarithmic sin- 
gularity for the specific heat C oc In |T — Td and critical 
indices /3 = 1/8 for the order parameter [m oc (Tc — T)^, 
T < Tc] and 7 = 7/4 for the staggered susceptibility 




-2 2 
(h - h^)/T 




FIG. 4: (Color online) Specific heat C for the frustrated bi- 
layer with Ji = 1, J2 = 5 {hi = 9 and /i2 =5): Exact diago- 
nalization data versus hard-square predictions. Upper panel: 
C in dependence on the hard-core parameter {h — hi)/T, cf. 
Eq. (|5.1I) [symbols: spin system with N ^ 20 {Af ^ 10), 
double-dashed line: hard-square model with Af — 10]. Lower 
panel: C in dependence on the temperature T for various 
values of magnetic field h [open symbols: spin system with 
A*' = 16 (A/" — 8), lines: hard-square model with Af = 8, black 
filled up-triangles: Monte-Carlo results for the hard-square 
model for large Af up to 800 x 800 for /i = 8.9]. 



[Xst oc |T — Tel"'']. This conclusion drawn from the clas- 
sical hard-square model taking into account only inde- 
pendent localized-magnon states (non-overlapping hard- 
square states) is reliable only for a quite small interval 
{hi — h)/hi <C 1. However, in the next section we will 
demonstrate that the range of validity of the classical 
description can be significantly extended including inter- 
acting localized-magnon states (overlapping hard-square 
states). 
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VI. LATTICE-GAS MODELS WITH FINITE 
REPULSION 

Taking into account also the interacting localized- 
magnon states, i.e., low- lying excitations, we will present 
a significantly improved lattice-gas description of the 
low-temperature thermodynamics of the considered spin 
models in this section. For that we use the information on 
the energies and degeneracies of these excitations, given 
in Sees. IIIII and llVI We start with the partition function 
for hard-core objects Zim(r, /i, N) and use its form given 
in Eq. (|5.2p . The hard-core rule is taken into account 
by the factor i?({nm}). To relax this rule we can pre- 
serve the form of the partition function, but we have to 
replace i?({n„}) by a modified factor g-^^/^^ ^f™") 
taking into account the excitation energy V ^ Ji. Then 
we arrive at the following formula instead of Eq. (|5.2p 



Z{T,h,N)^ Z^M{T,h,N) 



n 1—0,1 ?ijvr— 0,1 



(6.1) 



where ^ — hi — h and the sum X](mp) ''^i^^ over all 
nearest-neighbor bonds on the underlying lattice. The 
limit V/T — > oo now corresponds to the hard-core 
limit given in Eq. (j5.2p . since for V/T — >■ oo we get 
g-(WT)E(™p) _^ i?({n„J), i.e., the excitations get 
zero statistical weight. Note that in the improved lattice- 
gas description there is now an explicit temperature 
dependence in addition to the hard-core combination 
{hi-h)/T. Apart from the trivial factor e-'-'^™''''^/^^'^ 
the partition function Zlm(T, h, N) in (|6.ip is the grand- 
canonical partition function 5ig(T, ^,A/') of the lattice- 
gas model with finite nearest-neighbor repulsion < < 
00. Instead of Eq. (|5.3p we now have 



Fi^MiT,h,N) _ Efm h rinSig(T,Ai,AA) 



TV 



N 2 2 



with 



ni— 0,1 TiAT— 0,1 



where 



(6.2) 



(6.3) 



(6.4) 



in the one-dimensional case or 



A/'x K 

'H{{nm}) = E E [-/^""xm^ 



rrij-—! nil, — 1 



+v( )] (6.5) 

in the two-dimensional case. From Sig(r, fJ.,M) the ther- 
modynamic quantities can be found in usual way, cf. Ap- 
pendices |X] and |B] 



The following remarks are pertinent. Firstly, we no- 
tice that the initial quantum spin model (|2.ip has 2^ 
states and obviously not all of them are included in the 
effective models. Thus, the hard-core models contain 
either ip^^^ « 1.272^ states (one-dimensional case) or 
k(1)^/2 ~ 1.226^ states (two-dimensional case). The 
lattice-gas model with finite repulsion has 2^/^ w 1.414^ 
states. 

Secondly, we note that the particle-hole symmetry 
has some useful consequences. After making the trans- 
formation rim — ^ nm = 1 ^ nm in Eq. (|6.4p or 
rim^niy -> fim^rny = 1 ^ "-mxm« in Eq. wc arrive 

at the Hamiltonian H{{hm}) with + 2V instead of 
/i and shifted by Af{—fi + V) in the case (|6.4p or with 
—fi + W instead of /i and shifted by N{—^Ji -\- 2V) in 
the case (|6.5p . This fact implies, that Sig(T, /x, A/") = 
e^(''-^)/^Sig(r,-^+2F,A/') in the one-dimensional case 

or Sig(T,Ai,AA) = e^^^-^^^'^^i^{T,-^i + W,N) in the 
two-dimensional case. In particular, this yields identi- 
cal ground-state residual entropies at the fields hi and 
/i2. Moreover, the lattice-gas model with finite repulsion 
provides similar descriptions of the initial quantum spin 
model around both characteristic fields hi and /12. 

Thirdly, it is useful to introduce the on-site spin vari- 
ables (T = ±1 related to the site occupation numbers 
ri = 0, 1 as follows: cr = 2n — 1 and n = (1 + cr)/2. Then 
Eq. (|6.4p becomes the Hamiltonian of the antiferromag- 
netic Ising chain in a uniform magnetic field 



A/" 



m+1 ) 



rn—1 



F = ^i = hi-h, hi=2Ji + J2, J^^ (g.8) 

whereas Eq. (|6.5p becomes the Hamiltonian of the 
square-lattice antiferromagnetic Ising model in a uniform 
magnetic field 



n 



2 2 



-Fcr. 



m^—l niy — l 



F = ^-V, n^hi-h, hi^4Ji + J2, J =^ ^(7.) 

Let us mention again that for V ~ Ji we get correspon- 
dence to initial quantum spin systems. Note further, 
that the residual entropy present in the initial quan- 
tum spin systems at h = hi and /i = /12 corresponds 
to the known residual entropy of the Ising antiferromag- 
net at F = zk2J' (one-dimensional case) and F = zLAJ" 
(two-dimensional case).— From Ref. [33 we know that 
the ground state entropy per site at the critical fields 
is ln[(l -I- •\/5)/2] = 0.4812 . . . (one-dimensional case) and 
w 0.4075 (two-dimensional case) that coincides with the 
corresponding data for the considered quantum spin sys- 
tems reported in Sees. IV Al and lVBl 
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Now we discuss the low-temperature properties of the 
quantum spin models under consideration on the basis of 
Eqs. ()6.2p - (j6.7|) considering separately the frustrated 
two-leg ladder and the frustrated bilayer. 



A. Frustrated two-leg ladder 

The one-dimensional lattice-gas model with finite 
nearest-neighbor repulsion admits rigorous analysis. 
With the help of the transfer-matrix method we get 



-ze T ± 



\ 



ze T 



+ z 



(6.8) 



with z 



,{hi-h)/T 



Thermodynamic quantities for finite 



and infinite systems are given in Appendix |X1 Eqs. (|A11[) 
- (|XT4| and Eqs. (|M6)) - ([JT9|l . Obviously A± in Eq. 
(lelSj) transforms into A± in Eq. dO]) if V/T oo. It is 
interesting to note that according to Eq. (|6.8p A±(/i = 
2V) = e^/^A±(/i 0) and hence Sig(T,/^ = 2V,N') = 
e-'^^/^Sig(T, /i = OjA/"). This relation was mentioned 
already on the basis of particle-hole symmetry for hard- 
core objects. 

The transfer-matrix calculation in the case of different 
chemical potentials fiA and /is on the sublattices A and 
B leads to the following result for the grand-canonical 
partition function 



(6.9) 



[1 + ZA + ZB + ZAZse T 



+ZAZB 1 - e 



= 0. 



^± in Eq. (|6.9|) transforms into ^± in Eq. (|5.5p if 
V/T — > oo. With (|6.9p we can calculate the staggered 
susceptibility Xst{T,h,N), see Eqs. (jAlsp and (|A20j) in 
Appendix O 

We start with a general discussion of the low- 
temperature properties of the frustrated two-leg ladder 
based on its correspondence to the Ising chain (|6.6p . 
The one-dimensional Ising antiferromagnet exhibits an- 
tiferromagnetic long-range order along the line T = if 
|r| < 2 J and ferromagnetic long-range order along the 
line T = if |r| > 2J', whereas for any nonzero temper- 
ature it is in a disordered phase. In lattice-gas language 
this means that at T = the lattice is empty when /i < 0, 
one sublattice is completely occupied and the other one 
is empty (two-fold degenerate phase) when < n < 2V, 
and all lattice sites are occupied when 2V < fx. In terms 
of the initial quantum Heisenberg ladder this means that 
at T = the Ising-like antiferromagnetic long-range or- 
dered phase occurs if h2 = hi — 2Ji < h < hi only, 
i.e., for magnetic fields within the one-half magnetiza- 
tion plateau. Thermal fluctuations destroy perfect orders 



and a smooth crossover from the empty lattice to the lat- 
tice occupied by J\f localized magnons takes place at any 
fixed nonzero temperature as h decreases from above the 
saturation field hi to zero. 

Now we turn to numerics for finite systems. We fix 
Ji = 1 and set for concreteness J2 = 4 > J| « S.OOJi. 
Note that with increasing of J2 the lattice-gas descrip- 
tion is expected to become better, since excitations not 
described by the lattice-gas model are shifted to higher 
energies. In Fig. [5] we compare some results for the spe- 
cific heat of finite systems obtained from exact diago- 
nalization of the spin systems and from the lattice-gas 
formulas, see Eqs. ([Ell), (E^, dH]) and (KTT\ - (|M0)) 
with V = Ji = 1. For illustration we show in the upper 
panel also the hard-dimer result {V — >■ 00). We observe 
a very good agreement until T = 0.5, whereas the hard- 
dimer description is not appropriate at that temperature. 
For T = 1 deviations between exact diagonalization data 
and lattice-gas predictions become noticeable. However, 
increasing of J2 to J2 = 10 the exact diagonalization re- 
sult is again indistinguishable from the lattice-gas predic- 
tions (see dotted line and pentagons in the upper panel 
of Fig.©. 

It is important to note, that the obtained results for 
the thermodynamic quantities refer not only to finite sys- 
tems shown in Fig. [SJ but also to thermodynamically 
large systems. The thermodynamic quantities in the 
limit N 00 are given in the Appendix |X] by Eqs. (|A16[) 
- (|MOl) . Hence our findings for the thermodynamics of 
the frustrated two-leg ladder with J2 > together the 
ground-state analysis given Ref. |23 lead to a comprehen- 
sive description of that frustrated quantum spin model 
in the strong coupling regime. 



B. Frustrated bilayer 

Next we consider the lattice-gas model with finite re- 
pulsion, that is relevant for the frustrated bilayer (|2.ip . 
i.e., a lattice-gas of squares on the square lattice, where 
partial overlap is allowed, cf. Sec. IIVI For small fi- 
nite lattice-gas systems we use exact formulas for ther- 
modynamic quantities (see Appendix |B]). For large fi- 
nite lattice-gas systems we perform classical Monte Carlo 
simulations^ (see also Appendix IB|) . 

We start with a brief summary of the known results 
for the phase diagram of the corresponding square-lattice 
Ising antiferromagnet with nearest-neighbor exchange J7 
in a field F (|6.7|) ,— which sets the benchmarks in our 
further discussion. In contrast to the one-dimensional 
case, the two-dimensional model is known to have an an- 
tiferromagnetic long-range order within a restricted part 
of the half-plane "magnetic field F - temperature T" . A 
critical line separating the ordered regime along which 
thermodynamic quantities become singular has been dis- 
cussed in many papers.—^— Several closed-form formulas 
of the critical line T'c(r) were suggested and compared 
with numerical results. Clearly, along the line T = 
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FIG. 5: (Color online) Specific heat C for the frustrated two- 
leg ladder with Ji = 1, J2 = 4 {hi = 6, /12 = 4) and Ji — 1, 
J2 = 10 (/ii = 12, /i2 = 10): Exact diagonalization data for 
= 16 (TV = 8) versus lattice-gas predictions for TV = 8. 
Upper panel: C in dependence on the hard-core parameter 
{h~h\)/T for T — 0.5 and T = 1 [open symbols: spin system 
with J2 = 4, filled pentagons: spin system with J2 = 10, long- 
dashed and dotted lines: lattice-gas model, double-dashed 
line: hard-dimer model]. Lower panel: C in dependence on 
the temperature T for various values of magnetic field h [sym- 
bols: spin system with J2 = 4, lines: lattice-gas model] . Note 
that the short-dashed {h = 4.5) and the double-dashed line 
{h = 5.5) coincide because of the particle-hole symmetry in- 
herent in the lattice-gas model. For the spin model the cor- 
responding symbols (up-triangles and circles) also coincide at 
lower temperatures. 



the antiferromagnetic phase exists if |r| < AJ', whereas 
along the line F = the antiferromagnetic phase ex- 
ists below To/ J = 2/\n{V2 + 1) ^ 2.269 185 (Onsager's 
zero-field critical point^). For the lattice-gas model the 
corresponding critical line Tc{fJ,) is in the half-plane fj, - 
T, and we get Tq/V = l/[21n(v^-h 1)] « 0.567 296 at 
/i = 2V. The critical line crosses the /x-axis at /i = 



FIG. 6: (Color online) Phase diagram of the s = 1/2 Heisen- 
berg antiferromagnet in a magnetic field on the frustrated 
bilayer lattice with Ji = 1, J2 = 5 (i.e., hi = 9, /i2 = 5) in 
the half-plane "magnetic field h ~ temperature T" . The crit- 
ical line Tc{h) (solid line) separates the long-range ordered 
phase, in which the localized magnons occupy one sublattice 
and the other one is empty (below the curve) and the disor- 
dered phase with a random distribution of localized magnons 
(above the curve). The dashed line corresponds to the critical 
line according to the hard-square description. 



and /X = 4y . For the initial frustrated quantum Heisen- 
berg bilayer we have to set V = Ji. Then at T = the 
long-range ordered phase occurs if hi — 4Ji < h < hi, 
i.e., for magnetic fields within the one- half magnetiza- 
tion plateau h2 < h < hi. The corresponding critical 
line Tc{h) is in the half-plane h - T, and we get the max- 
imal critical temperature Tq = Tc[h = (hi — h2)/2] — 
Ji/[21n(y2 + 1)] w 0.567 296Ji. 

Hence, in the frustrated bilayer we have various pos- 
sibilities to pass from the disordered to the long-range 
ordered phase (where localized magnons occupied only 
one of the two sublattices): (i) Fixing the field h, /12 < 
h < hi, and decreasing of temperature T to T < Tc(h). 
(ii) Fixing the temperature T,T < Ji/[2 ln(y2-|-l)], and 
decreasing of h starting from above hi. (iii) Fixing the 
temperature T, T < Ji/[2ln{V2 + 1)], and increasing 
of h starting from below /i2. Crossing the critical line 
Tc{h) the critical behavior is that of the two-dimensional 
Ising model. In Fig. [5] wc show the phase diagram of the 
frustrated bilayer which is a retranslation of the corre- 
sponding phase diagram of the square-lattice Ising anti- 
ferromagnet in a field .■^"ia Note, however, that we have 
reproduced this phase diagram by our Monte-Carlo sim- 
ulation of the classical lattice-gas model (for more details 
see below). 

A short remark about the hard-square case, which de- 
scribes relevant physics of the frustrated bilayer (|2.ip 
around hi and small T, is expedient here. The criti- 
cal line as it follows from the hard-square model reads: 
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Tc(/i) = {hi - h)/\B.Zc with Inz. « 1.3340 (the dashed 
hne in Fig. H]), see Refs. [Hill The critical behav- 
ior which emerges while crossing the curve Tc{h) — 
{hi — h)/ \nzc for the hard-square model also belongs to 
the two-dimensional Ising model universality classJ^ii^ 
The results shown in Fig. |6] demonstrate that the hard- 
square phase diagram coincides with the lattice-gas phase 
diagram only around the point h ~ hi and T < 0.3. 

To estimate the validity of the lattice-gas phase dia- 
gram for the quantum Heisenberg antiferromagnet on the 
frustrated bilayer lattice ()2.ip we again compare exact 
diagonalization data for finite bilayer spin systems with 
the lattice-gas predictions, see Fig. [71 It is also useful to 
compare these results with corresponding ones obtained 
using the hard-square model, see Fig. |4l From the re- 
sults reported in Fig. [7] one concludes that the lattice-gas 
model provides a very good description of the considered 
finite quantum spin system with Ji = 1, J2 = 5 in a 
wide range of magnetic fields at least up to temperatures 
about T = 0.5 Moreover, if J2 acquires a large value, 
J2 = 10; lattice-gas predictions remain very good even 
at T = 1.0 (see the upper panel in Fig. [7]). Hence, we 
have evidence that the phase phase diagram presented in 
Fig. [6] is indeed valid for the s = 1/2 Heisenberg anti- 
ferromagnet in a magnetic field on the frustrated bilayer 
lattice. 

Let us discuss two further aspects of the data shown 
in Fig. [T] (i) Comparing lattice-gas results (long-dashed 
line in the upper panel) with hard-square results (double- 
dashed line in the upper panel) the limited tempera- 
ture range of validity of the hard-square picture is ob- 
vious, (ii) The extra low-temperature maximum in C{T) 
is present in a wide range of magnetic fields. It is well 
described by the lattice-gas model. 

After having demonstrated the quality of the lattice- 
gas description of the low-temperature thermodynamics 
for small systems we consider now thermodynamically 
large systems. We have used classical Monte Carlo sim- 
ulations for the lattice-gas model with finite repulsion 
which reproduce reliably the low-temperature proper- 
ties of the frustrated bilayer in a wide range of mag- 
netic fields. Note, however, that for the special value 
of /i = {hi + /i2)/2 = 2Ji + J2 we face the zero-field 
square-lattice Ising model and hence in this limit we 
have a set of analytical equations for thermodynamic 
quantities known from Onsager's solutionJ^i^ In Fig. [5] 
we show temperature dependences of the specific heat 
C, staggered susceptibility Xst, and the entropy S ob- 
tained from classical Monte Carlo simulations.— Due to 
the particle-hole symmetry (which we have confirmed ex- 
plicitely by our Monte-Carlo calculations) the tempera- 
ture dependence is identical at fields h = h2 + Ah and 
h = hi — Ah. As a main feature, clearly seen for large sys- 
tems for h2 < h < hi, a. divergence of the specific heat C 
(Fig. m upper panel) and the staggered susceptibility Xst 
(Fig. ISl middle panel) appears at a critical temperature 
where the order-disorder phase transition takes place, cf. 
the phase diagram shown in Fig. [51 As already discussed. 
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FIG. 7: (Color online) Specific heat C for the frustrated bi- 
layer with Ji = 1, J2 = 5 {h\ = 9, h2 — 5) and Ji — 1, 
J2 = 10 {hi = 14, /i2 = 10): Exact diagonalization data ver- 
sus lattice-gas predictions. Upper panel: C in dependence on 
the hard-core parameter {h - hi)/T for N ^ 20 (Af = 10) 
and T = 0.5 and T = 1 [open symbols: spin system with 
J2 — 5, filled pentagons: spin system with J2 — 10, long- 
dashed and dotted hues: lattice-gas model, double-dashed 
line: hard-square model]. Lower panel: C in dependence on 
the temperature T for various values of magnetic field h for 
TV = 16 {Af — 8) [symbols: spin system with J2 = 5, lines: 
lattice-gas model]. The lines for h — 6 and h — 8 coincide be- 
cause of the particle-hole symmetry of the lattice-gas model. 



the critical behavior belongs to the two-dimensional Ising 
model universality class, see Sees. IV Bl and IVIBI Note 
that the low-temperature maximum in C(T) found for 
smaller systems (see Fig. [71 lower panel) is masked by 
the logarithmic divergence. However, for magnetic fields 
near /12 and hi a characteristic bump below the diver- 
gence occurs. The high-temperature maximum in C(T) 
present for the spin system cannot be described correctly 
by the lattice-gas model. Another feature worth to be 
mentioned is the behavior of the staggered susceptibility 
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FIG. 8: (Color online) Specific heat (upper panel), staggered 
susceptibility (middle panel), and entropy (lower panel) ver- 
sus temperature for the frustrated bilayer with Ji = 1 and 
J2 = 5 {hi = 9 and /12 = 5) for different magnetic fields ob- 
tained from Monte Carlo simulations for the effective lattice- 
gas model with finite repulsion V = Ji. 



Xst for T 0. While there is an exponential decay of 
Xst to zero for h2 < h < hi, precisely at h = hi and 
/i = /i2 it diverges as 1/T. With respect to the tempera- 
ture dependence of the entropy shown in the lower panel 
of Fig. [8] it is worthwhile to make the following remark. 
To obtain the entropy S(T, h, N) we perform integration 
according to Eq. (|B8p . Note, however, that in Eq. (jB8|) 
the contribution at T = 0, i.e., S{T — 0,h,N), is not 
included. This is correct if h ^ hi oy h ^ h2- How- 
ever, for h = hi OT h = h2 there is a nonzero ground- 
state residual entropy which is therefore missed in the 
corresponding curve (diamonds) in the lower panel of 
Fig. [8l Taking into account the constant of integration 
S{T ^ 0,hi,N) ^ S{T ^ 0,h2,N) = ln[«;(l)]iV/2 leads 
to a shift of the curves S{T, hi,N)/J^ and S{T, /i2, N)/J^ 
upward by sa 0.4075 and recovers a monotonic depen- 
dence of temperature profiles S{T) as h varies from 4.0 
to 5.5, see Fig. [51 lower panel. However, if the magnetic 
field is close to hi or /i2 (see crosses for h = 4.9 and pen- 
tagons for h = 4.99 in Fig. [5]) as a remnant of the residual 
entropy present for h = hi and /i = /12 the entropy re- 
mains large up to quite low temperatures T ~ |/i — /12I or 
\h — hi\, respectively. 



VII. CONCLUSIONS 

In the present paper we have demonstrated that the 
thermodynamic quantities of two particular quantum 
many-body systems, namely the frustrated ladder and 
bilayer s = 1/2 Heiscnberg antiferromagnets in a mag- 
netic field h (|2.1[) . can be obtained via classical lattice- 
gas-model calculations in a wide range of magnetic fields. 
For the one-dimensional ladder model by means of the 
transfer-matrix method even a complete analytical de- 
scription is possible, whereas for the two-dimensional bi- 
layer model well elaborated classical Monte Carlo sim- 
ulations can be used. The reason for this significant 
simplification lies in the simple structure of low-energy 
levels of the quantum spin system which emerges due to 
frustrations in the strong-coupling regime. The ground- 
state magnetization curve in this regime exhibits plateaus 
at zero magnetization and at one-half of the saturation 
magnetization. The classical lattice-gas model leads to 
an excellent description of the quantum spin models up 
to quite large temperatures of the order of the exchange 
constants in the field region of the one- half plateau, i.e., 
for h2 < h < hi, a,s well as magnetic fields slightly below 
/i2 and above hi. 

Some prominent features of the considered systems are 
as follows: a ground-state residual entropy at h ~ h2 and 
h = hi that may be of particular interest for magnetic 
cooling^"— liSiiiiiid^ and a second order phase transition 
at a critical temperature Tc{h) > 0, h2 < h < hi found 
for the two-dimensional bilayer system, where the critical 
behavior is that of the two-dimensional Ising antiferro- 
magnet. 

Finally, we mention that lattice-gas approach elabo- 
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rated in the present paper can be extended to similar 
(although different) models, such as the frustrated three- 
leg ladder— in one dimension or the bilaycr systems con- 
sisting of two triangular or honeycomb lattices. We leave 
the discussion of these models for further studies. 
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Appendix A: One-dimensional lattice-gas models 

In this appendix we collect some analytical results for 
the one-dimensional (i) hard-dimer model and (ii) lattice- 
gas model with finite nearest-neighbor repulsion obtained 
by means of the transfer-matrix method. These formulas 
can be used to calculate the relevant physical quantities 
at low temperatures for the s = 1/2 frustrated Heisen- 
berg two-leg ladder. 

We start with the one-dimensional hard-dimer model, 
see Eqs. (|5.4p and (|5.5p . For the entropy, the specific 
heat, the average number of hard dimers, the uniform 
susceptibility, and the staggered susceptibility after sim- 
ple but lengthy calculations we get 



SiTV^N) _ I f ^ (Inz)z X+^-'-X^^-' 



(Al) 



+ 



(In z)2 



Af 



(AA-1): 



VT+Tz 



7V-1 \ A/"-! 



C{T,h,N) _ (Inz)^z / _ 2z \ A+^~' - A 



1 + 42/ A+^+A 



A+^ + A_-^ 



(A2) 



M{T,h,N) _ n n _ 



z A+-^-i - A_^-i 



N 



Af' Af VT+Tz A+-^ + A_^ 



(A3) 



Tx{T,h,N) 



VTTT 



(AA-i): 



VT+Tz 

M-2 ^ y^_U-2 



1 - 



2z 



l + 4z 



X+-^-' - A_-^-i 
A,^+A_^ 



A+^ + A_^ 



-Af\ 



A+^ + A_^ 



(A4) 



Txst{T,h,N) _ 2 A+^-i + A_^-i 



Af 



2 A+-~ + A 



, \ A/' 



z X+-^-' - X^^-' 
2VT+T X+^ + A_^ ' 



(A5) 



respectively. We recall that here z ~ eS'^^ '»)/^^ Inz = 
{hi - h)/T, and A± = (1 ± Vl + 4z)/2 [see Eq. (lO)) ]. 

In the limit A/" — > oo the formulas (|Aip - (jA5[) become 
much simpler 



SiT,h,N) 



In A, 



(Inz)z 1 

VTT4iA7' 



(A6) 



C{T,h,N) _ (lnz)2z 



Af 



(l + 4z)^ 



A/(T, h,N) _^ n n 



Af Af' Af VTTTX+'' 

TxiT,h,N) 



Af 



(H-4z)^ 



(A7) 

(A8) 
(A9) 
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TXstjT, h, N) _ z (AlO) specific heat, the average number of hard dimers, the 

M Vl + 4z uniform susceptibihty, and the staggered susceptibihty 

after simple but lengthy calculations we get 

We turn to the one-dimensional lattice-gas model with 
finite repulsion, see Eqs. (|6.8p and (|6.9p . For the entropy. 



S(T,KN) 1,/, AA^.^X 2(lnz)z - (lnz^)u;(l - w) X+^-' X^^-' 1 , A+^l^l+A-^ , ^ ^ . 

- :ln A+ +A_ = , A/- , , A, -i\nw)w ^ ^ \ ^ (AH) 



Af Af \ ^ J 2^(1 -wy+Az A+-^ + A_^ 2^ ' A+^ + A_ 

CiT,h,N) b+A+-^-^+b_A_-^-^ ^^ a^^A+^-^ + «^_A_-^-^ ^ g+A^-^-^ + a_A_-^-^ V 

" A+-^ + A_-^ ^ A+-^ + A_-^ -^1, A+-^ + A_^ )' 

2(lnz)2 — (lnK7)w(l — u') 1 

2V(l-.0^+4z 2(1""^"' 

^.-■^(^-)^--(^^-)^-(^-^-)t[^(^-)-(^"-)-(^:-)]' I (A12) 
2V(l-u02+4z 2[(l~i«)2 + 4z]^ 2^ 

M{T,h,N) _ n_ n_ 2z ~ w{l ~ w) A+-^'^ - A_-^'^ w A+-^'^ + A_-^~^ 

M "^"AA' 2^/O^WTTz A+^ + A_^ ^2" A+^ + A_^ ' ^ 

rx(T,^Ar) _ rf+A+^-^+rf_A_-^-^ , ^^ 4A+^-^ + c^_A_^-^ ^r ( ^l^^llf^L^Y 

A+^ + A_-^ +(-^-1) XT^TT:^ -^I^ A+^ + A_^ ) ' 

2z — 1^(1 ^ w) w 



c± =± 



2^(1 -w)2 +4z 2 



l2 



2z — w(l — 2^) [2z — wfl — w)!^ w , . , 

d± = ± — ^ , — =F ^ + — , (A14) 

2^(l-i(;)2 + 4z 2[(l-«;)2-h4z]^ 2 

rxst(r,fe,jv) _ z A+^-^ + A_-^-^ z A+-^-^ - A_-^-^ 

2(1 + A+^ + A_-^ 2^/(1 -ti')2+4z A+^ + A_^ ' ^ ' 



respectively. We recall that here z = e^''^ '»)/'^^ Inz = 
{hi - h)/T, X± = [l + w± A/(l-w)2-h4z]/2 [see Eq. 
(|6.8p ] and we have also introduced the notations w = 
g(fti-/i-,/i)/T^ Inw ^(hi~ h - Ji) /T. Evidently in the 
limit w Q Eqs. (jAlip - (jAlSp transform into Eqs. 
dHI) - jM]). 

In the limit A/" oo the formulas (|All|) - (jAlsp be- 
come much simpler 



S{T,h,N) 



lnA4 



2(lnz)z — (lnu')z«(l — It;) 1 



2^(l-u;)2-H4z 

C{T,h,N) _ 
N " A I 



-(In w)?!! 



A2^ 



^,(A16) 



(A17) 



n 



M{T, h, N) 
AA 

2z — w{l — w) 
2^/(1 -i«)2 +42 



A/" A+ Al 



1 - 


n 




w 


1 


2" 


K 



Txst{T,h,N) 



(1 -I- w)yj{l - wY -I- 4z 



(A18) 



(A19) 



(A20) 



Again after inserting -u; = into Eqs. (|A16p - (jA20[) we 
obtain Eqs. (|X6]) - (|A10|) as it should be. 
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Appendix B: Two-dimensional lattice-gas models 

In this appendix we collect some formulas for the two- 
dimensional (i) hard-square model and (ii) lattice-gas 
model with finite nearest-neighbor repulsion which we 
use in our direct calculations of thermodynamic quanti- 
ties for small systems and classical Monte Carlo simula- 
tions for large systems. 

We can obtain thermodynamic quantities for finite 
hard-square models by direct calculations starting from 
the definition of the grand-canonical partition function 

K 

2 

Shc(^,AA) = JZhc(n,AA)z" (Bl) 

and knowing the canonical partition functions Zhcin,Af) 
for n = 0, 1, . . . , A/'/2 (calculation of these numbers are 
feasible for small J\f). Really, Shc(-Z, A/") is a polynomial 
of order M/2 and calculations of thermodynamic quan- 
tities are doable although rather tedious. Thus, for the 
entropy, the specific heat, the average number of hard 
squares, and the uniform susceptibility we find 

S{T,h,N) ^ ln[Shc(2,A/')] - (Inz)n, (B2) 
C(T, h, N) = (In (Ti^-TT^y (B3) 



M{T,h,N) = --n, (B4) 

Tx{T,h,N)=^-if, (B5) 

where z = e'''^^''^/-^ is the activity and (...) = 
E^io ^hc(»^,A/')2:"(. . .)]/Shc(2;,A/') denotes the grand- 
canonical average for the considered finite hard-square 
model of N sites. 

For completeness we give here canonical partition func- 
tions for some finite lattices: Zhc(".,8) = 1,8, 12,8,2 for 
n = 0,1,2,3,4, Zhc(n, 10) = 1,10,25,20,10,2 for n = 
0, 1, 2, 3, 4, 5, Zhc(n, 16) = 1, 16, 88, 208, 228, 128, 56, 16, 2 
for 71 = 0,1,2,3,4,5,6,7,8. 

To calculate the staggered magnetization 

M,t{T,h,N)^rrA-nE (B6) 
and the staggered susceptibility 

~2{nAnB - riA ub) (B7) 

we have to introduce the quantities Zhc(".A7 "-s; A/") 
which are the numbers of spatial configurations of n 
hard squares, where ha of them occupy the sublat- 
tice A and ns = n — ua of them occupy the sublat- 
tice B. Obviously, Zhc(ri,A/') = X)Ii^=o ^hc("-A, A/"). 



We have to refine the definition of the grand-canonical 
average making it sensitive to the sublatticc indices. 
When the staggered component of the acti vity van- 
ishes, i.e., za = zb = z, we have (...) = 
[En=o ^"E:,=o^hc(n^,nB;AA)(. . .)]/Shc(z,AA). 

We give here (n^ , ; A/") for some finite lattices: 
Zhe(l,0;8) = 4; Zhc(2,0;8) = 6, Zhe(l,l;8) = 0; 
Zhc(3,0;8) = 4, Zhc(2,l;8) = 0; Zhe(4,0;8) = 1, 
Zhc(3,l;8) = Zhc(2,2;8) = for AA = 8, Zhc(l,0;10) = 
5; Zhe(2,0; 10) = 10, Zhe(l, 1; 10) = 5; Zhe(3,0; 10) = 10, 
Zhc(2,l;10) = 0; Zhc(4,0;10) = 5, Zhc(3,l;10) = 
Zhe(2,2;10) = 0, Zhc(5,0;10) = 1, Zhc(4,l;10) = 
Zhc(3,2;10) = for A/'= 10. 

Formulas (|B3[) - (|B7p are also used for obtaining 
Monte Carlo predictions for large hard-square systems. 
In Monte Carlo simulations we calculate rTJ, tTb, n^, 
n^, and riAnB for a given z = e^^^"'^^^"^. As a result, 
we obtain n = rTX + TTb and hence M{T, h, N) (jB4[) and 

— = 71^ — tTa^ + n\ — TTb^ + 2{riAriB ~ ^nAns) and 
hence x{T,h,N) (|B5]) and C{T,h,N) (|B3|. Then the 
entropy is obtained by integration 

S{T, h, N) = £ dT' ^^^''^'^\ (B8) 

Moreover, Monte Carlo data yield \M^t{T,h,N)\ 
and Xst{T,h,N) 

We turn to the two-dimensional lattice with finite 
nearest-neighbor repulsion. Starting from the formula for 
the grand-canonical partition function Sig(r, A/") ()6.3p 
and the definition of the grand-canonical average 

^lg(J ,//,A/) 

we immediately get 



5(T,/7,iV)=lnSig(r,A*,AA) + ^^^^i^, (BIO) 

C(r, /7, N) = ^^({""'})^^^({"-})! (Bll) 

for the entropy and the specific heat, respectively. For the 
uniform magnetization, the uniform susceptibility, the 
staggered magnetization, and the staggered susceptibility 
we formally have the same expressions as in Eqs. (jB4p . 
(jB5|) . (jB6[) . and (|B7|) . however, with the grand-canonical 
average defined in Eq. (|B9|) . In the limit V/T — >■ oo 
Eqs. (|BT0)) . ((BTTI) transform into Eqs. (HI]), dBS]) since 
Eq. (|B9p becomes the grand-canonical average for the 
hard-square model and 'H{{nm}) ~^ —fJ-n, where n is the 
grand-canonical average number of hard squares. 

We use formulas (|BIO)) . (|BTT1) . ((B4)) . (|B^ . ([Be]) , and 
(jB7|) with the grand-canonical average (|B9p for direct 
calculations for small finite systems encoding easily the 
required computations for A/" = 8, 10 in a short Fortran 
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program. F or large sy ste ms wc obtai n from Monte Carlo 
simulations H{{nm}) and 'H'^{{n,m}) yielding the specific 
heat C{T,h,N) (jBll[) and by integration the entropy 
S{T, h,N), see Eq. (jB8[) . Furthermore, we also compute 



nA, TiBi n^, and n^n^ to obtain the magnetizations 
and the susceptibilities, see Eqs. (|B4|) . (|B5|) . (|B6p . and 
(iBTl) . 
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